#include <kuai/mol/molio.hpp>
#include <kuai/sim/trajec.hpp>

#ifndef _KUAI_SIM_GROMACS_HPP_
#define _KUAI_SIM_GROMACS_HPP_

namespace kuai {

	class GromacsTOPReader 
		: public SpaceReader
	{
	public:
		virtual bool update(std::istream& stream, PhaseSpace& space) const;
	};


	class GromacsTrajecPlayer
		: public TrajecPlayer
	{ 
	public:
		DECLARE_WORKER_CREATOR(GromacsTrajecPlayer)

		GromacsTrajecPlayer();
		virtual size_t current();
		virtual void update(PhaseSpace& space, Index tid=0, Index nThreads=1);
		virtual void open();
		virtual void close();

	private:
		struct GromacsTRNHeader /* This struct describes the order and the	*/
			/* sizes of the structs in a trjfile, sizes are given in bytes.	*/
		{
			bool	bDouble;    /* Double precision?                            */
			int		ir_size;	/* Backward compatibility		        */
			int		e_size;		/* Backward compatibility		        */
			int		box_size;	/* Non zero if a box is present			*/
			int		vir_size;   /* Backward compatibility		        */
			int		pres_size;  /* Backward compatibility		        */
			int		top_size;	/* Backward compatibility		        */
			int		sym_size;	/* Backward compatibility		        */
			int		x_size;		/* Non zero if coordinates are present		*/
			int		v_size;		/* Non zero if velocities are present		*/
			int		f_size;		/* Non zero if forces are present		*/

			int		natoms;		/* The total number of atoms			*/
			int		step;		/* Current step number				*/
			int		nre;		/* Backward compatibility		        */
			RealNumber t;			/* Current time					*/
			RealNumber lambda;		/* Current value of lambda			*/
		} ;

		struct t_energy {
			RealNumber	e;		/* The current energy.					*/
			double		eav;    /* The running average       				*/
			double		esum;	/* The sum of energies until now.			*/
			RealNumber	e2sum;	/* The sum of the square of energies until now		*/
		} ;


		struct t_enxframe {
			RealNumber t;	            /* Timestamp of this frame		              */
			int      step;          /* MD step				              */
			int      nre;           /* Number of energies			      */
			int      ndisre;        /* Number of distance restraints	              */
			int      nblock;        /* Number of following energy blocks              */
			std::vector<int> nr;           /* Number of things in additional blocks (nblock) */
			int      e_size;        /* Size (in bytes) of energies		      */
			int      d_size;        /* Size (in bytes) of disre blocks	              */
			int      nr_alloc;      /* Allocated size of nr and block                 */
			int      e_alloc;       /* Allocated size (in elements) of ener           */
			int      d_alloc;       /* Allocated size (in elements) of rav and rt     */
			int      *b_alloc;      /* Allocated size (in elements) of each block     */
			t_energy   *ener;         /* The energies                                   */
			RealNumber *disre_rm3tav; /* Time averaged data for distance restraints     */
			RealNumber *disre_rt;     /* Instantaneous data for distance restraints     */
			RealNumber **block;       /* Additional energy blocks (nblock x b_alloc[b]) */
		} ;

	private:
		bool nextTRN();
		bool nextEDR();
		bool parseEDRNames();

		bool readInteger(FILE* f, int& result) {
			bool ok = fread(&result, sizeof(int), 1, f) == 1;
			if (ok && _reversed) {
				result = swapbytes(result);
			}
			return ok;
		};
		bool readString(FILE* f, String& result) {
			int n = 0;
			if (readInteger(f, n)) {
				result.resize(n);
				return fread(&(result[0]), 1, n, f) == n;
			}
			return false;
		};
		bool readFloat(FILE* f, float& result) {
			bool ok = fread(&result, sizeof(float), 1, f) == 1;
			if (ok && _reversed) {
				result = swapbytes(result);
			}
			return ok;
		};
		bool readDouble(FILE* f, double& result) {
			bool ok = fread(&result, sizeof(double), 1, f) == 1;
			if (ok && _reversed) {
				result = swapbytes(result);
			}
			return ok;
		};
		bool readFloat(FILE* f, float result[], size_t n) {
			bool ok = fread(result, sizeof(float), n, f) == n;
			if (ok && _reversed) {
				for (size_t i = 0; i < n; ++i) {
					result[i] = swapbytes(result[i]);
				}
			}
			return ok;
		};
		bool readDouble(FILE* f, double result[], size_t n) {
			bool ok = fread(result, sizeof(double), n, f) == n;
			if (ok && _reversed) {
				for (size_t i = 0; i < n; ++i) {
					result[i] = swapbytes(result[i]);
				}
			}
			return ok;
		}

		FILE* file;
		GromacsTRNHeader _trnHeader;
		bool _reversed;
	};

}

#endif
